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Abstract 

Wavelet  bases  are  used  to  generate  spaces  of  approximation  for  the  resolution 
of  bidimensional  elliptic  and  parabolic  problems.  Under  some  specific  hypothe¬ 
ses  relating  the  properties  of  the  wavelets  to  the  order  of  the  involved  operators, 
it  is  shown  that  an  approximate  solution  can  be  built.  This  approximation  is 
then  stable  and  converges  towards  the  exact  solution.  It  is  designed  such  that 
fast  algorithms  involving  biorthogonal  multi  resolution  analyses  can  be  used  to 
resolve  the  corresponding  numerical  problems. 

Detailed  algorithms  are  provided  as  well  as  the  results  of  numerical  tests  on 
partial  differential  equations  defined  on  the  bidimensional  torus. 
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Introduction 


Variational  approximation  methods  for  partial  differential  equations  are  based 
on  weak  formulations  and  on  suitable  spaces  of  approximation. 

Wavelets  are  known  to  be  unconditional  bases  for  a  large  variety  of  spaces 
and  therefore  are  good  candidates  for  the  generation  of  approximation  spaces 
for  partial  differential  equation  problems.  The  goal  of  this  paper  is  to  show  that 
moreover,  wavelet  bases  may  lead  to  fast  and  adaptive  numerical  resolution  of 
the  corresponding  approximations. 

In  this  paper,  as  in  previous  papers  (J.  Liandrat  and  P.  Tchamitchian  [10], 
[11]),  the  wavelets  are  used  to  expand  the  approximated  solution  of  a  partial 
differential  equation  as  well  as  to  approximate  the  kernel  of  the  differential 
operator.  They  are  not  used  only  to  perform  “the  linear  algebra”  (G.  Beylkin 
[2])  related  to  more  classical  methods  of  resolution. 

Starting  with  an  expansion  of  the  form  /  =  Ylx  <  /?  the  solution 

of  the  equation  Lu  f  where  L  is  a  constant  coefficient  elliptic  differential 
operator  is  =  f  f(p)K(x,z/)d^  where  K(x,y)  = 

Under  suitable  conditions  that  will  be  made  precise  later,  the  functions 
as  well  as  are  pseudo  wavelets,  very  close  to  wavelets  (Y.  Meyer 
[13]).  This  turns  out  to  provide  a  stable  approximation  of  u.  However,  the  effi¬ 
ciency  of  the  corresponding  numerical  approximation  of  u  relies,  at  least  in  this 
work,  on  the  hierarchic  structure  of  multiresolution  analysis  since  it  provides 
fast  tree  algorithms.  We  will  show  that,  if  the  operator  satisfies  suitable  condi¬ 
tions  that  will  be  made  explicit  later,  then  the  above  mentioned  pseudo  wavelets 
are  directly  related  to  biorthogonal  multiresolution  and  wavelets.  Under  these 
conditions,  competitive  numerical  algorithms  involving  0(N)  or  0(N  log  N)  op¬ 
erations  can  then  be  derived. 

This  paper  provides  the  analysis  of  the  problem  and  exhibits  the  correspond¬ 
ing  numerical  schemes.  It  is  then  organized  as  follows. 

The  first  part  is  devoted  to  the  general  concept  of  biorthogonal  multiresolu¬ 
tion  analysis  on  In  the  second  part  we  focus  on  the  problem  of  the  sta¬ 

bility  of  the  multiresolution  framework  under  the  action  of  constant  coeflScient 
elliptic  operators.  The  cases  of  homogeneous  and  non  homogeneous  operators 
are  treated  separately.  The  third  part  deals  with  the  numerical  algorithms  while 
the  last  section  is  devoted  to  numerical  tests  related  to  the  resolution  of  elliptic 
and  parabolic  equations  in  bidimension al  spaces. 

I  Generalities:  Biorthogonal  Multi-Resolution 
Analysis  in  L^{TR^) 

The  concept  of  multiresolution  is  at  the  basis  of  our  construction  and  we  there¬ 
fore  start  with  a  short  description  of  it  : 
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Definition  I.l  (Y.  Meyer,  [12]) 

A  r-regular  multiresolution  analysis  of  L^{1R^)  is  a  sequence  of  increasing 
closed  subspaces  Vj,  j  €  Z,  Vj  C  satisfying  the  following  conditions 

=  is  dense  ini  W; 

ii)  f{x)  E  Vj  f{2x)  E  Vj+i; 

Hi)  f{x)  E  Vo  f{x  —  k)  E  Vq,  VAr  E 

iv)  there  exists  a  function  $  in  Vq,  such  that  the  set  of  functions  {^{x-k),  k  E 

is  a  Riesz  basis  ^  for  Vo; 

v)  the  function  $  is  regular  and  localized  :  $  is  is  almost 

everywhere  differentiable,  and  for  almost  every  x  E  for  every  integer 
a  <r  and  for  all  integer  p,  ii  exists  Cp  such  that 

\d^^x)\<Cp{l^\x\)-P  .  (1) 

A  consequence  of  ii),  Hi)  and  iv)  is  that  each  Vj  is  generated  by  the  family 
of  functions  ~  k),  k  E 

For  simplicity  reasons,  we  will  only  consider  the  case  n  =  2,  but  all  the 
results  presented  to  this  article  can  be  generalized  in  any  dimension.  We  will 
always  use  for  vectors  a  contracted  notation:  if  e  is  a  bidimensional  vector  then 
e  =  (ei,e2). 


I.l  Orthogonal  multiresolution  analysis 

To  build  an  orthonormal  multiresolution  analysis,  the  Riesz  basis  {^(.  —  k),kE 
is  first  orthonormalized  in  such  a  way  that  the  resulting  orthonormal  basis 
is  still  of  the  form  {$(.  —  k),kE  Z^]. 

The  wavelets  are  introduced  via  the  orthogonal  complement  of  Vo  in  Vi :  Wo- 
More  precisely,  if  E  is  the  set  of  all  vertices  in  the  unit  cube  [0, 1]^,  and  if 
E*  =  E  \  {0},  we  have  the  following  theorem: 

Theorem  LI  There  are  3  functions  e  E  E*,  in  Wo,  such  that  the  collection 
{^^{x  —  k),  k  E  Z"^ ,e  E  E""]  is  an  orthonormal  basis  ofWo.  Moreover,  each 
satisfies  the  same  property  (1)  of  regularity  and  localization  as  $  and,  moreover, 
satisfies  the  following  cancelation  property 

3  m  E  JIN  ,  such  that  V  Ar  =  (fei,  ^2)  €  0  <  <  m, 

f  . 

/  =  0  . 

Jm? 

collection  of  vectors  {e^,  A  6  A},  in  a  Hilbert  space  H  is  a.  Riesz  basis  if  any  vector 
X  £  H  can  be  written  in  a  unique  way  as  a:  =  ^  where  is  finite  and 

defines  a  norm  equivalent  to  ||a:||jj. 
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The  scaling  invariance  property  ii)  implies  that,  for  all  j,  the  family k  E 
e  G  £'*}  is  an  orthonormal  basis  of  Wj.  We  will  also  use  the  following 
contracted  notation:  A  G  A^}  where  Aj  =  {A  =  2”^  4*  |)  ,  A?  G  G 

Indeed,  there  is  a  straightforward  bijection  between  Kj  and  the  set  of  pairs 
{(£,  k),  k  G  G  £'*}.  We  will  also  use  the  following  sets:  A  =  Uj^^Aj  and 

A  n  »  jm  A  . 

From  the  inclusion  Vb  C  Vi  the  following  scaling  relation  can  be  derived: 

$(x)  =  -  0  . 

lez^ 

while  from  Wo  C  Vi  one  obtains  the  following  detail  relation: 

r  (ar)  =  -  0>  ^eeE*  . 

It  is  very  useful  to  transform  these  relations  using  the  Fourier  transform 
which  is  given  by  the  equality 

m  =  I  f{x)e-^^^dx. 

JlR 


Indeed,  the  scaling  relations  then  become 


$(0  =  Mo(^/2)%/2) 

(3) 

and 

r(0  =  M,(^/2)$(^/2) 

(4) 

where  Mo(0  =  E/ez"  and  M,{^)  =  are  (7°° 

odic  functions. 

This  leads  to  : 

27r  peri- 

oo  .  f  f 

$(o = 

j=i  j=i 

(5) 

The  following  conditions  are  satisfied 

^  +  7re)Me-(^  +  :re)  =  6,,,  ,  V  (e,e')  € 

(6) 

and 

Mj(7re)  =  V(£:,e)  €  E^ 

(7) 

and  are  called,  following  electrical  engineering  terms,  the  quadrature  mirror  filter 
conditions.  We  will  also  call  the  functions  M.  quadrature  mirror  filters. 
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Conversely,  it  is  shown  in  A.  Cohen  et  ah  [4]  that  four  27r  periodic  functions 
s  £  E  satisfying  the  quadrature  mirror  filter  conditions  (6)  generate 
through  (5)  an  orthogonal  multiresolution  analysis  if  some  specific  conditions 
are  satisfied. 

Remark: 

♦  In  this  paper,  we  will  often  use  specific  multi-resolution  analyses  of  L^{IR^) 
based  on  a  tensorial  product  of  multiresolution  analyses  of  L^{1R).  More 
precisely,  such  analyses  are  defined  as  follows:  if  (Vj)  is  the  sequence  of 
spaces  of  a  ID  multiresolution  analysis  and  if  Wjj<p,  mo  and  mi  are  re¬ 
spectively  the  related  wavelet  spaces,  the  scaling  function,  the  associated 
wavelet  and  the  quadrature  mirror  filters,  then,  the  sequence  of  spaces 
(h^  ),  defined  as  Vj  =  Vj  ®  Vj  is  a  multiresolution  analysis  in  IR^.  More¬ 
over,  =  (p{xi)(f{x2)  is  the  corresponding  scaling  function;  (Wj),  with 
Wj  =  wavelet  spaces;  the  three  generating  wavelets  are 

and  =  V’(a^i)^(^2); 

Me(0  mei(^i)me2(^2)  with  e  £  E  ene  the  quadrature  mirror  filters. 

1.2  Biorthogonal  Approach 

A  relaxation  of  some  properties  of  orthogonal  multiresolution  analysis  can  be 
performed  using  the  one  biorthogonal  approach.  This  approach  provides  some 
flexibility  since  it  allows  to  distribute  the  relevant  properties  of  the  multires¬ 
olution  (number  of  zero  moments,  compact  support  or  regularity)  to  the  two 
involved  multiresolution  analyses.  Moreover,  it  will  turn  out  to  be  that  the 
biorthogonal  framework  is  “stable”  under  the  action  of  a  large  class  of  opera¬ 
tors  while  the  orthogonal  framework  is  “fragile” . 

Definition  L2  We  call  biorthogonal  muliiresolution  analysis  of  two 

multire solution  analysis  (Uj)j^z  such  that  there  exists  two  families 

of  corresponding  scaling  functions  r  and  r  such  that:  {vjk,  Tjfk')  =  bjjfSkk'  for 
oil  j,  f  G  and,  k  and  k'  £  . 

In  this  case  we  define  the  wavelets  spaces  Xj  and  Xj  as  Uj+i  =  Uj  ^  Xj, 
=  Uj  0  Xj  with  UjEXj,  Ujl.Xj,  and  we  introduce  the  functions  = 
2*^ .  —  k)  and  ^a  =  2-^ ^^(2-^  .  —  k),  e  £  E,  that  generate  respectively  Xj  and 
Xj  and  such  that  (^a?  =  S^efbjjtSkk'^ 

Moreover,  following  the  construction  of  orthogonal  multiresolution  analysis 
we  define  2x4  filters  (i.e,  (7°^  27r  periodic  functions),  and  e  £  E, 
associated  with  the  two  biorthogonal  multiresolution  analyses.  These  filters 
satisfy  the  biorthogonal  quadrature  mirror  filter  relations  equivalent  to  (6)  that 
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are 


V  5  and  e'  E  E,  (  E  [0,  27r]^  , 

P,{^  +  ire)Pe’{^  +  Tre)  =  6es’  ,  (8) 

e^E  _ 

Pq^^TTC  )  =  ej^£2,e2  • 

As  in  the  orthonormal  case,  the  generalizations  of  relations  (3)  and  (4)  relate 
the  scaling  functions  and  the  wavelets  to  the  filters  as 

oo 

m  =  Pomm^)  =  =  Psimmni  (9) 

i=i 

oo  ^ 

?(o  =  PoimHm  =  ^  (o  =  Pcmfim-  (lo) 

i=i 

As  in  section  Ll,Jhe  question  to  know  is  under  which  conditions  the  biorthog- 
onal  filters  and  satisfying  the  quadrature  mirror  filter  conditions  (8)  de¬ 
fine  two  biorthogonal  multiresolution  analyses?  Again,  a  specific  condition  is 
required  and  has  been  formulated  in  A.  Cohen  et  aL  [4].  We  will  use  a  weaker 
version  of  this  formulation  adapted  to  the  case  of  functions  with  fast  decay.  It 
can  be  expressed  in  term  of  the  following  theorem  for  which  a  complete  proof 
can  be  found  in  Pj.  Ponenti  [16]  and  W.  Dahmen  and  A.  Michelli  [5]. 

Theorem  L2  Lti(x>  Q  and  lei  Pe{^)  and  Pe(0?  e  E  E,  be  eight  27r  periodic 
functions  satisfying  the  biorthogonal  quadrature  mirror  filters  conditions  (8), 
Defining  r,  f,  6^  and  9  using  formula  (9)  and  (10),  if 


-  there  exist  C  and  /i  >  0  such  that  for  all  ^  E  IR^ 


|r(OI  <  C(l  +  |e|)-l-^ 

i?(oi  <  c-a + lei)-'-'' , 

(11) 

-^k  ,etk'eZ\ 

{t{x  -  k),  t(x 

II 

1 

(12) 

-  and  if 

L  ifp  ^  ' 

L  i«p 

(13) 

then 

♦  the  sequences  of  subspaces  {Uj)j£Z}  <^'^d  (P})j€Z  generated  respectively 
by  {vjkjk  E  and  {rjk^k  E  are  two  biorthogonal  multiresolution 
analyses; 


♦  the  wavelet  families  {  9\{x)  =  —  k),  A  G  A  }  and  {^a(3^)  = 

2^9^{2^x  —  k),  A  G  A}  are  two  biorthogonal  Riesz  bases  of  L^{1R^), 
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II  Constant  Coefficient  Elliptic  Differential  Op¬ 
erators  and  Biorthogonal  Multiresolution  Anal¬ 
ysis 

II.  1  General  results 

The  starting  point  of  this  section  is  the  following  remark.  Given  ($a)  a  family 
of  orthonormal  wavelets  and  knowing  /  =  the  solution  of  the 

equation 

Lu  =  f  ,  (14) 

where  L  is  an  elliptic  operator  of  order  s  is,  at  least  formally, 

u  =  5](/,^a)L-M^a]  .  (15) 

When  L  and  L~^  are  bounded  on  the  families  and  {L*$a} 

are  two  biorthogonal  Riesz  bases. 

The  question  we  address  now  is  related  to  what  happens  in  the  specific  case 
of  wavelets  when  the  operator  L  is  unbounded  (as  in  the  case  of  a  differential 
operator) . 

In  the  following  paragraphs,  we  first  see  that,  assuming  some  compatibil¬ 
ity  conditions  between  and  L,  the  two  families  of  functions  and 

are  wavelets  or  pseudo  wavelets  (Y.  Meyer  [13])  depending  on  whether 
the  operator  is  homogeneous  or  not.  Then,  we  show  that  in  some  cases,  a 
biorthogonal  framework  embedding  and  {T*^a}  can  be  built. 

To  be  more  precise,  we  take  (Vj)  an  r-regular  multiresolution  of  L^{IR?)  con¬ 
structed  using  a  tensorial  product  of  two  ID  multiresolutions  and  L  a  constant 
coefficient  elliptic  differential  operator.  Let  us  write  L  =  Z)o<a<5  where 

D  is  the  operator  We  define  in  a  standard  way  the  symbol  of  L  as 


It  is  a  polynomial  in  ^  of  degree  s  and,  if  /  €  C°°(iR^),  we  have  the  well  known 
formula 


(17) 

We  note  formally 

0^  and^"  = 

(18) 

and  more  generally. 

e{  =  and  91  = 

(19) 

Note  that  when  L  is  homogeneous,  ^^(ar)  =  2^6"{2^x 
2^9^{2^ X  —  k)  while,  in  general  this  is  not  true. 

Then  we  have 

—  k)  and  0\{x)  = 
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Theorem  ILl  Given  a  family  of  r  regular  wavelets  with  m  +  1  zero  moments, 
if  L  is  a  homogeneous  operator  or  if  L  is  an  inhomogeneous  operator  with  a 
strictly  positive  symbol  of  order j  >  0,  and  such  that  r  >  s  and  m>  s,  then 
Regularity:  ^  G  and  0^  E 

Localization:  for  all  multi-indices^  7'  and  7  such  that  |7'|  <  r  —  s  and 
I7I  <  r  -f-  5  —  1  and  all  integers  I  E  IN ^ 

[  \^'6\{x)\  <  Cy  2JlT''l2J'  (l  +  2Jla;-A|)"' 

I  (l  +  2J>  - 


Cancelation:  Let  x* 


—  {Xi\  X2^), 

X^e^{x)dx  =: 
x^0^{x)dx  = 


0, 

0, 


0  <  |^r|  <  m  +  a 
0  <  \k\  <  m  —  a  , 


where  a  =:  s  in  the  homogeneous  case  and  a  =  0  in  the  inhomogeneous  one. 


A  complete  proof  of  this  theorem  can  be  found  in  Y.  Meyer  [13]  and  Ph. 
Tchamitchian  [17]. 

Remark: 


♦  Following  Y.  Meyer  [13]  and  according  to  Theorem  ILl,  a  factorization 
of  the  operators  L*  and  L~^  can  be  performed  as  L*  =  C  o  and 
zz  C  oT^  where  F  is  a  diagonal  operator  in  the  wavelet  basis  defined 
as:  rpx  ^  2-^Va  and  where  C  and  C  are  bounded  on  Lp"  and  defined  by 


C  :  Ox.  C  :  ^x-^  Ox.  (20) 

The  operators  C  and  C  act  just  as  a  transformation  between  two  bases. 
The  operator  F  is  nothing  else  but  the  classical  preconditioning  operator 
for  elliptic  problems  (S.  Jaffard  [8])  that  mimics  a  diagonal  derivation  in 
the  wavelet  basis. 

In  other  words,  thanks  to  this  factorization,  the  computation  of  the  image 
of  a  function  by  an  elliptic  operator  or  its  inverse  can  be  transformed  into 
a  well  conditioned  problem  using  a  diagonal  operator  in  a  suitable  wavelet 
basis.  This  is  essential  since  it  provides  the  numerical  stability  of  the 
further  developed  algorithm. 


We  can  then  rewrite  (15)  as 


^We  call  a  multi-index  any  couple  of  integer  7'  =  (7i,72)and  d^' 0  = 


(21) 


An  important  issue,  as  far  as  numerical  applications  are  concerned,  is  the 
computation  of  the  sum  (21).  Indeed,  even  if  this  sum  corresponds  to  a  pseudo 
wavelet  decomposition,  fast  algorithms  for  the  computation  of  the  sum  are  not 
available.  In  the  framework  of  this  paper,  the  fast  algorithms  are  linked  to 
the  concept  of  multiresolution  analysis  presented  in  section  I.  We  prove  in 
the  following  sections  that,  under  suitable  condUions,  the  construction  of  a 
multiresolution  analysis  embedding  the  function  0\  is  possible.  Moreover,  we 
provide  explicit  expressions  of  the  quadrature  mirror  filters  required  for  the  fast 
implementation  of  (21). 

The  starting  point  of  our  construction  is  due  to  P.G.  Lemarie  [9]  who  con¬ 
structed  in  the  one  dimensional  case  two  biorthogonal  multiresolution  analyses 
from  an  original  orthogonal  one  and  from  the  derivata  operator.  We  generalize 
this  approach  to  any  dimensions  and  for  any  homogeneous  elliptic  operator. 

In  contrast  to  the  classical  constructions  of  multiresolution  analysis,  this  is 
an  inverse  problem.  Indeed,  knowing  the  two  dual  wavelet  bases  we  can  define 
two  sequences  of  subspaces  {Xj)  and  (Xj): 

Xj  ^  span'^^x )  ^  j  ^  ^  1  f 221 

Xj  =  span{^A,  A  e  A/,  /  <  j}  , 

The  open  question  is  the  following:  how  can  we  construct  two  sequences  of  sub¬ 
spaces  (Uj)  and  (Uj)  for  which  (Xj)  and  (Xj)  play  the  role  of  two  biorthogonal 
wavelet  spaces? 

In  other  words,  the  problem  is  Jhe  construction  of  the  generalized  scaling 
functions  r  and  r  related  to  0  and  0. 

In  the  case  of  non  homogeneous  elliptic  operators  the  approach  used  in  the 
homogeneous  case  can  not  be  transposed  and  we  will  not  be  able  to  define  a 
multiresolution  framework  embedding  the  space  sequences  (Uj)  and  (Uj)-  How¬ 
ever,  we  will  show  that  the  essential  property  of  embedding  spaces  as  well  as 
the  existence  of  scaling  (3)  and  detail  (4)  relations  can  be  saved.  This  will  allow 
us  to  derive  fast  and  stable  algorithms  to  sum  up  (21)  even  in  the  case  of  non 
homogeneous  operators. 

11.2  The  case  of  homogeneous  elliptic  operator 

In  that  case  the  natural  candidates  for  r:  Q'l'e  not  defined  in  for  the 

basic  reason  that  ^  does  not  belong  to  the  range  of  T,  or,  in  other  words,  suffers 
from  a  lack  of  zero  moments.  Using  a  preconditioning  operator,  we  will  adapt 
the  function  $  to  the  operator  (i.e.  we  will  transform  $  so  that  it  enters 
the  range  of  I)  while  preserving  the  two  scale  relations  (3)  and  (4).  We  will 
then  check  that  the  resulting  multiresolution  analysis  is  fitted  to  the  functions 
6x  and  Ox  previously  defined. 

More  precisely,  we  have  the  following  theorem: 
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Theorem  IL2  Lei  (Vj)  be  the  family  of  embedded  spaces  of  an  r -regular  mul- 
iiresoluiion  analysis  ofL^{lR^),  lei  L  be  a  homogeneous  operator  of  order  s  and 
of  symbol  a  and  lei  S  be  a  2tt  periodic  function,  not  vanishing  on  ]0,27rp  and 
equivalent  to  a  in  zero.  If  m  >  s+l,  r  >  s  +  1,  and  if  the  eight  functions  Pe{^) 
and  ^  ^  E  defined  as 

PoiO  =  2^  ^  ^  ^ 

P.(0  =  2*^  M,(0,  RiO  =  (24) 

are  C^,  a  >  0  then  they  satisfy  the  quadrature  mirror  filter  conditions  (8)  and 
define  two  biorihogonal  multiresoluiion  analyses.  The  corresponding  wavelets 
are  the  functions  6  and  0  and  the  scaling  functions  r  and  r  are  derived  following 
(9)  and  (10). 

Proof: 

□  By  construction  all  the  filters  are  with  some  a  >  0  and  they  satisfy 
the  biorthogonal  quadrature  mirror  filter  conditions  (8). 

The  only  point  to  prove  is  the  convergence  in  L^{1R^)  of  the  infinite  prod¬ 
ucts  (9)  and  (10)  defining  the  two  scaling  functions.  We  will  use  the  following 
lemma: 

Lemma  ILl  Let  p{x),  x  E  IB?,  be  a  homogeneous  polynomial  of  degree  s, 
and  let  S  and  C,  be  27r  periodic  functions;  then  the  following  propositions  are 
equivalent: 


ii) 


nc(2->x) 

j=i 


p{x)  ’ 


(25) 


Six)  =  2^Cix/2)Six/2)  , 
5(a;)  ~  p(a;)  . 

r— ►U 


(26) 

(27) 


Proof: 

□  The  equation  (26)  is  obtained  from  (25)  written  for  x  and  x/2,  while  (27) 
is  derived  from  (25)  when  x  0  since  necessarily  (7(0)  =  1. 

Conversely,  (25)  is  obtained  from  (26)  .  Indeed,  since 


Sjx) 

p{x) 


Cix/2) 


S{x/2) 

pix/2) 


S{2-^x) 

p{2~^x) 


N 


nc(2-^'^) 

i=i 
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thanks  to  (27)  we  obtain  (25)  when  N  ^  oo.  m 

This  lemma  allows  us  to  calculate  the  infinite  product  (9)  and  (10),  and  to 

get: 


CO  CO  -  oo 

j=l 

OO  00  OO 

%)  =  n  ^0(2“''^)  =  n<^(2-^On”^o(2-^^i)mo(2-^'6)  , 

;=1  j  =  l  i  =  l 

and  finally, 

^  "  W) 

The  function  $  being  r-regular  the  conditions  (11)  are  trivially  satisfied. 
Defining  the  function  ^  et  ^  by  (9)  and  (10)  we  get 

^^(0  =  P.(?/2)r(^/2)  = 

f(0  =  Am?(^/2)  =  |||. 

The  wavelet  admissible  condition  (13)  is  immediately  satisfied  thanks  to  Theo¬ 
rem  ILl.  Finally  the  assumption  (12)  is  satisfied  by  construction,  that  completes 
the  proof.  ■ 

Remarks: 

•  Note  that  here,  thanks  to  homogeneity,  the  subscript  a  recovers  its  classical 

“wavelet”  meaning  since  in  that  case,  ^a(^)  =  x-k)  said  9\{x)  = 

2^^^e{{2^x-k). 

•  The  relation  (25)  is  a  generalization  in  any  dimension,  of  the  classical 
formula 

ncos(2-^x)  =  ^ 

J=1 

used  by  P.G.  Lemarie  for  the  first  order  differential  operator  in  one  di¬ 
mension. 

•  The  function  S  can  be  interpreted  as  the  symbol  of  a  difference  operator 
that  we  will  call  Dl.  If  5  is  a  trigonometric  polynomial,  then  Dl  is  a 
finite-difference  operator  and  S  is  .  The  fact  that  S{x)  (t{x)  is 

just  the  translation  that  Dl  is  consistent  with  X.  Moreover,  5  removes 
exactly  the  singularity  of  r  for  =  0.  Conversely  for  r,  the  singularity 
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given  by  S  at  points  (27rn,27rn),  n  £  will  be  removed  exactly  by  a  in 
zero  and  by  the  zeros  of  ^  at  points  (27rn,  27rn),  n  £  ZS* ,  Notice  that  this 
last  point  won’t  be  true  for  inhomogeneous  operators. 

From  a  certain  point  of  view,  Dl  can  be  seen  as  a  preconditioner  for  L 
since  r  and  r  are  defined  by 

T  =  L-^Dl^  =  (30) 

•  In  one  dimension  there  is  a  canonical  choice  for  S  and  therefore  for  Di 
such  that,  if  the  function  $  and  ^  have  a  compact  support,  then  r,  r, 
6,  and  9  are  also  compactly  supported.  Indeed,  in  that  case  we  have 
necessarily  cr(^)  =  a  €  (T.  Therefore,  the  canonical  choice  for  S 
is  5(0  =  ^i~^Y  and  Dl  is  then  a  non-centered  finite- 

difference  approximation  of  L  of  order  1.  Indeed,  it  is  well  known  (see 
I.  Daubechies  [6])  that  the  quadrature  mirror  filters  related  to  orthogonal 
compactly  supported  functions  ^  and  ’F  are:  mo(0  =  (l  +  >C(0 

and  mi(0  =  (l  —  £(^  +  tt)  where  £  is  a  finite  trigonometric 

polynomial.  Then  we  get 

PoiO  =  r  (1  +  m,  _ 

PliO  =  (l  -  (l  -  £(4  +  tt), 

PoiO  =  2-^  (1  +  (1  +  -^(0, 

PiiO  =  2-^a-U^e^^  (1  -  - 

which  proves  that  Pq,  Pi,  Pq  and  Pi  are  also  finite  trigonometric  polyno¬ 
mials.  Then,  using  the  following  lemma  borrowed  from  G.  Deslauriers  and 
S.  Dubuc  [7]  we  deduce  that  the  functions  r,  r,  9,  and  9  have  compact 
support. 

Lemma  II.2 

um  =  En^;v,7ne--«  wHh  =  I,  nr=ir(2-i0  «  an 

entire  function  of  exponential  type.  In  particular,  it  is  the  Fourier  trans¬ 
form  of  a  distribution  with  support  in  [A^i,  N2]^ 

Clearly,  this  canonical  form  is  no  longer  available  if  the  space  dimension 
is  larger  than  1  since  the  multidimensional  quadrature  mirror  filters  can 
not  be  factorized  as  above. 

II.3  The  case  of  inhomogeneous  elliptic  operator 

Here  the  non  homogeneous  property  of  the  operator  is  obviously  not  adapted 
to  the  scale  invariance  property  of  the  multiresolution  analysis.  We  will  see 
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however,  that  introducing  at  each  level  a  new  scaling  function,  an  embedded 
family  of  spaces  can  be  constructed  which  preserves  the  mathematical  properties 
relevant  for  numerical  applications. 

The  natural  candidates  for  t\  are  a  the  functions  They  are  well 

defined  in  L’^{1E?)  but  suffer  now  from  a  lack  of  localization  when  j  increases. 
Indeed,  we  have 

lim  II  [^.^]  (a;)  -  G{x  -  2"^'^)  ||  =  0 

j— ^+oo 
— ^xo 

where  G  is  the  Greens  function  of  the  operator  L  defined  as  G{^)  =  For 

example,  when  L  =  1- A,  G(^)  =  1/(1+^^)  and  G{x)  =  G  decreases  fast 

but  mathematical  and  numerical  difficulties  come  from  the  fact  that  the  family 
of  functions  {G{x  —  k2~^),  k  G  is  not  a  set  of  functions  rescaled  with  j  (in 
other  words,  this  family  is  not  obtained  by  rescaling  and  translation  of  a  single 
initial  function).  This  implies  that  the  control  of  the  localization  by  the  index 
j  is  lost.  It  follows  that  the  family  k  E  is  not  a  good  basis  to 

reconstruct  our  solution. 

Let  us  show  now  that,  however,  a  process  very  close  to  the  one  used  in 
the  homogeneous  case  will  provide  an  efficient  algorithm  for  the  summation  of 
formula  (21). 

We  mimic  the  construction  performed  in  the  homogeneous  case.  Let  (Vj) 
be  an  r-regular  multiresolution  analysis,  let  L  be  an  elliptic  operator  of  order  s 
with  constant  coefficients  and  a  its  symbol  (we  now  suppose  cr(^)  >  ctq  >  0  V(J), 
Let  us  also  define  the  homogeneous  polynomial  of  order  s,  (j,  as  the  principal 
part  of  <7,  and  let  S{^)  be  a  27r  periodic  function  with  5(<f)  ^  ^  ^  where 

n  will  be  fixed  later. 

Then,  Vj  E  we  define  a  difference  operator  Dj  by  its  symbol  S{^/2^). 

Following  the  previous  section,  we  define  Vj,  k  £  Z 

Tjk  =  (31) 

and 


(32) 


By  construction,  and  thanks  to  the  fact  that  L  is  a  constant  coefficient 
operator  we  have  Tjk{x)  =  Tj{x  —  ^)  and  6j}^{x)  =  dj{x  —  where 


and 


(33) 


(34) 
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Remark 


♦  The  functions  tj  mimic  the  function  f  defined  in  (28).  Unfortunately,  it  is 
not  possible  to  define  the  equivalent  of  r  (28)  since  L^jk  i  Note 
however  that,  by  chance,  (15)  involves  directly  the  9\  functions. 

Then,  with  Pq  and  P,  defined  in  (23)  and  24  we  get  the  following  scaling 
and  detail  relations: 

liO 

Remark: 

•  An  important  point  is  that  the  filters  Pe  are  independent  of  the  scale 

index  j  as  it  is  originally  the  case  for  standard  multiresolution  analysis. 
Furthermore  ,  since  they  are  defined  by  (23)  and  24,  the  filters  P^  are 
directly  related  to  the  homogeneous  operator  L  of  symbol  <7  if  n  =  s. 
This  point  is  essential  since  it  means  that  if  Do  is  consistent  to  T,  then 
the  tree  algorithms  related  to  the  multiresolution  spaces  (Uj)  and  used  to 
sum  up  (21)  are  stable  even  if  the  functions  tj  (x)  are  not  standard  scaling 
functions.  Indeed,  5(^)/cr(^)  0  ^nd  therefore  tj(0)  =  0. 

In  other  words,  even  if  the  functions  Tjk{x)  are  used  as  scaling  functions  on 
the  range  of  L,  they  have  zero  moments  as  wavelets  have. 

Finally,  we  can  prove  the  following  theorem: 

Theorem  IL3  For  0<n<s,  s<r  and  s  <  m,  the  functions  Tjk  defined 
by  (33)  and  (31)  satisfy 

\d'^rjk{x)\  <  C;  (1  +  2^>  -  .  (36) 

Ifn^s  and  ifPe{^),  £  E  E  are  then  stable  tree  algorithms  are  available 

30  <  C  <  C"  <  +00  such  that  if  f  —  dx6x,  then 

C^\dx\^  <\\ff<  c>  (37) 

A  A 


=  2Pom^+^)fj+i{0 
=  2n(^/2^'+')?+i(0 .  Veei;*. 


Proof: 

•  Since  5  €  the  we  can  apply  Theorem  11. 1  that  proves  the  localization 
inequality  (36). 
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•  The  relations  (35)  defines  directly  the  tree  algorithms  required  to  compute 
the  coefficients  {cjk}  such  that 


and  since,  when  n  =  s,  the  involved  quadrature  mirror  filters  could  be 
related  to  the  homogeneous  operator  T,  the  stability  of  the  algorithm  is 
equivalent  to  the  stability  of  the  transform  {cjk}  /  for  /  G  0^=0 


•  This  transform  is  stable  if  and  only  if  the  family  {rjk,  Ar  €  is  a  Riesz 
basis  of  0£Zo 
We  have 


II /IP 


\smnf  X 


Since  is  a  Riesz  basis,  since  5'(^)  is  bounded,  and  since  a  is  bounded 

from  below,  then  ||  /  |p<  C'Y^\cjk\^,  which  is  the  second  part  of  the 
inequality  (37). 

To  prove  the  first  part  of  (37)  we  use  again  the  fact  that  the  filters  Pq  and 
Ps  are  related  to  L. 

Indeed,  if  we  define  9\  replacing  a  by  <t  in  (29)  and  if  we  define  /  as 
/  =  dx6\  then  the  transform  f  {dx}  is  stable.  Moreover  thanks  to 

theorem  II.l,  the  operator  /  /  which  can  be  also  defined  as  ^x  ^ 

is  also  bounded  (Y.  Meyer  [13]).  Therefore  the  operator  /  {cIa}  is 
bounded  that  is  the  first  part  of  (37)  and  that  completes  the  proof. 


Ill  Approximation  and  Numerical  Resolution 
of  Elliptic  Problem  on  the  Torus 

This  section  is  devoted  to  the  approximation  of  elliptic  problems  on  a  sequence 
of  embedded  Galerkin  spaces  associated  with  a  multiresolution  analysis  and  to 
the  corresponding  numerical  algorithms. 
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Classically,  we  will  consider  the  problem  with  periodic  boundary  conditions 
to  avoid  the  difficulties  of  general  boundary  conditions.  We  will  use  a  r-regular 
multiresolution  analysis  of  the  torus  T[o,i]2  =  {R/2ZY  constructed  using  a  clas¬ 
sical  periodization  technique  (Y.  Meyer  [12]).  We  take  as  granted  that,  with 
minor  modifications,  the  results  proved  on  the  whole  line  can  be  transposed  to 
the  torus.  In  that  section,  homogeneous  and  inhomogeneous  operators  will  be 
treated  similarly. 


III.l  General  formulation 

The  general  formulation  of  the  problem  is 


Find  u  G  ^[o,i]2  suck  that 


i'P)  { 


Lu  =  f  (38) 

with  f  G  ^[o,i]2  dTid  L  a  constant  coefficient  elliptic  operator  of 
order  s. 


Standard  variational  approximation  (P.A.  Raviart  and  J.M.  Thomas  [14]) 
leads  us  to  look  for  the  solution  of  a  weak  problem  in  so  called  Galerkin  approx¬ 
imation  spaces  14  7  where  e  is  a,  scale  related  to  14  with  T[o,i]2. 

A  natural  choice  for  I4  is  I4  =  ^  where  belongs  to  a  multiresolution 

analysis  of  T[o,i]2  of  the  type  described  above.  Indeed,  we  then  have  the  following 
inequality  guaranteed  if  the  involved  multiresolution  analysis  is  r-regular, 

V  s  <  r,  3  c>  0  V/  G  <  C  2'^^  ||/||//.  (39) 


where  Hyv,  j  <  0,  stands  for  the  orthogonal  projection  on  Vf. 
Then  a  standard  Galerkin  approximation  writes 


Find  Up  G  C  ^[0,1]=  that 


(V)  { 


UyvLIlyvUp  =  Hyvf. 


(40) 


where  Hyp,  j  <  0,  stands  for  the  extension  operator  from  to  T[o^i]2. 

This  approach  leads  us  to  replace  L  by  the  approximation  UyvLIiyv.  The 
corresponding  numerical  algorithms  are  reduced  to  linear  system  solvers  once  a 
basis  of  Vp  has  been  chosen  (P.A.  Raviart  et  al.  [14]). 

§The  symbol  stands  for  the  periodization  operator  on  [0,  1]^.  We  recall  that  dim  14^  = 
2^^  =  1/3  dim  Wj  and  that  =  spanf^oo  =  !}• 
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Numerically,  an  optimal  choice  for  the  expansion  basis  is  the  wavelet  basis 

since  the  corresponding  stiffness  matrix,  is  sparse 
and  can  be  easily  (i.e  using  diagonal  matrices)  uniformally  preconditioned  (S. 
Jaffard  [8]). 

III.  2  A  different  approximation 

The  main  purpose  of  this  paper  is  to  define  a  different  approximation  of  L  and 
the  corresponding  numerical  algorithm.  Taking  the  set  of  functions  A  G 
(defined  in  (31))  and  defining  =  spanlL’^^oo^  ^  ^ 
an  approximation  of  u  as 

U,=  Y1 •  (41) 

AeAS 

Indeed, 

Up  =  L-^Uvv  f  =  f  (42) 

where  is  the  projection  on  orthogonal  to  . 

This  formulation  of  Up  shows  that  the  convergence  of  Up  towards  u  when 
p  — oo  is  straightforward  since  the  set  is  a  family  of  Galerkin  spaces  for  the 
suitable  space  of  definition  of  u. 

Moreover,  the  stability  of  the  algorithm  is  a  direct  consequence  from  the 
classical  preconditioning  properties  of  wavelet  base  expansions  (S.  Jaffard  [8]). 

IIL3  General  scheme 

The  numerical  algorithm  derived  from  the  previous  section  is  now  presented 
in  its  collocation  version.  We  call  Ij  the  set  of  points  (Z  f][0,2^[)^.  Then, 
Jj  =  k  e  Ij}  is  the  two  dimensional  regular  grid  of  scale  2~^  related  to 

[o,ip. 

For  the  numerical  implementation,  we  assume  that  the  space  Vp  is  such  that 
any  continuous  function  f  is  unambiguously  defined  by  its  values  on  the  set 
of  points  Jp.  This  assumption  (satisfied  by  the  even  order  spline  multiresolution 
we  will  use  in  the  numerical  tests)  allows  us  to  define  the  collocation  projection, 
Cyv,  from  the  set  of  real  sequences  {fk)kelp  as 

'^{fk)k^ip,Cvv{{fk)keij,)  =  f  <^f^Vp  and /(2“^Ar)  = /fe,  e  Ip. 

As  soon  as  we  define  f  E  its  coordinates  on  the  basis  the  opera¬ 

tor  Cyv  appears  as  a  discrete  convolution  operator  involving  a  so-called  inter- 
polant  filter  l$p(0*  The  operator  Cyr  is  also  a  convolution  operator  called  the 

^Again,  periodized  framework,  Aj  =  {A  =  2~^  (/:+  |)  ,/:  € 

Zf][0,V[)^,eeE*}, 
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point  value  operator  and  involves  the  point  value  filter  defined  from 

Obviously,  stnd  P14p(^)  are  inverse.  Let  us  remark  however, 

that  the  point  value  operator  always  exists  as  soon  as  the  functions  are 
continuous. 

For  the  implementation,  we  therefore  replace  by  Cyv  in  (42)  and  define 
therefore  Up  as 

Up  ^  L  ^Cy'pf  • 

Given  the  point  values  of  /  on  the  grid  points  Jp  ,  the  algorithm  provides  the 
values  of  Up  on  the  same  grid.  More  precisely,  the  algorithm  can  be  presented 
as  follows: 

1.  The  input  of  the  procedure  is  the  set  of  values  {f{Jp))  from  which  the 
interpolant  function  fp  E  is  constructed  using 

fp  =z  ^  Cpk^pk 

kelp 

2.  fp  is  then  decomposed  into  the  wavelet  subspaces  <  j  <p  — I  and 

Ko^as 

/=  E  (/’^a)^a+coo 

where  cqo  =  Uyfif)- 

3.  Up  then  becomes 

AGAr' 

where  6^  =  Va)^-  Here,  Cqq  =  coo/(t{0)  for  nomhomogeneous 

operators.  For  homogeneous  operator  Cqq  should  be  given.  Note  that  in 
that  case  cr(0)  =  0  and  /  should  have  at  least  s  vanishing  moments;  the 
fact  that  c'oo  should  be  given  corresponds  to  the  ill  posed  property  of  the 
initial  problem  in  L^. 

4.  Up  is  then  expanded  in  terms  of  the  set  of  functions  {^ph^  k  ^  Ip)  using 
the  tree  algorithms  related  to  and  ^pk 

Up  =  ^  ^pk'^pk  ‘ 

kelp 

5.  Finally,  the  grid  point  values  of  Up  on  Jp  are  estimated  using  the  point 
value  filter  PV- 
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It  appears  clearly  that  various  precalculations  should  be  performed.  In  the 
first  step,  the  interpolant  filter  related  to  ^pk  must  be  known;  for  the  second 
step,  the  orthogonal  multiresolution  analysis  quadrature  mirror  filters  should 
be  used  and,  for  the  fourt^tep  the  corresponding  biorthogonal  multiresolution 
quadrature  mirror  filters  are  required;  finally,  the  point  value  filter  related 
to  Tp  is  used  for  the  last  step. 

To  be  more  precise,  we  have  to  make  some  remarks  that  help  to  reduce  the 
complexity  and  storage.  For  steps  one  and  two,  tensorial  properties  can  be  used 
in  a  very  classical  way  to  reduce  the  2D- algorithms  involved  in  2(dim(l^^))^/^  x 
ID-algorithms.  It  is  then  enough  to  know  the  one-dimensional  interpolant  filter 
related  to  and  the  one  dimensional  quadrature  mirror  filters  ^  6  E. 

For  steps  four  and  five,  where  the  filters  and  the  point  value  filter  PV-r  are 
involved,  we  can  note  apply  this  simplification  and  we  have  a  full  2D-problem. 
We  are  now  able  to  summarize  all  these  precalculations  in  the  following  step 

0: 


0.  The  computation  of  the  following  filters  is  performed  (this  is  presented  for 
the  spline  multiresolution  analysis  case): 

-Interpolant  filter  related  to  (ppk->  I(pp'  analytical  formulas  in  one  di¬ 
mension  are  available  in  V.  Perrier  and  C.  Basdevant  [15]. 

-Orthogonal  ID  multiresolution  analysis  filters  m^:  analytical  formu¬ 
las  are  also  available. 

-Filters  Pe  of  the  Biorthogonal  Multi  resolution  Analysis:  These  filters 
are  constructed  from  and  formula  (23)  and  (24).  In  fact  only  P(o,i) 
and  P(i,i)  have  to  be  computed  since  we  have  P(i,o)(^ij^2)  =  P(o,i)(C25<^i)- 

-Point  value  filter,  PV~  related  to  This  filter  is  computed  from 

Tpk 

formula  (28)  and  the  analytical  expression  of  5,  a  and  4(^).  We  have 
successively 

T^{x)  =  l/(27r)  ®  ,  x  €.  B?  , 


?;(x„)  =  1/(2^)  Y,  E  +  2"'-)  €  Jp  .  (43) 

Practically,  43  is  truncated  according  to  a  prescribed  precision.  This  is 
possible  because  7>(0  decreases  fast. 

Remarks: 

•  One  should  again  emphasize  that  the  entire  algorithm  is  based  on  con¬ 
volution  operators.  Thanks  to  the  periodic  boundary  conditions,  the 
convolutions  can  either  be  performed  directly  or  using  a  discrete  Fourier 
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transform.  The  implementation  presented  in  this  paper  uses  the  Fourier 
transforms  since  it  is  optimal  for  non  compactly  supported  filters  on  non 
adapted  spaces  of  approximation. 

III.4  Detailed  Algorithm 

This  section  is  devoted  to  the  structure  of  the  code.  Basic  tools,  such  as  Fast 
Fourier  Transforms,  Convolution/ Decimation  algorithms,  or  Term  by  term  mul¬ 
tiplications  are  not  described. 

As  can  be  seen  from  the  general  scheme  presented  above,  the  main  code 
involves  only  two  more  elaborate  routines  that  will  be  called  the  Precalculus 
routine  (step  0),  and  the  tree  algorithm  routines.  The  tree  algorithm  routines 
may  or  may  not  use  the  tensorial  structure.  They  will  be  called  consequently  2D 
Tensorial  Tree  Algoriihm-D  (steps  2)  and  2D  Non  Tensorial  Tree  Algoriihm-I 
(step  4)  where  -D  and  -I  stand  for  direct  and  inverse.  We  recall  that  the  steps 
1  and  5  are  convolutions  and  the  renormalization  performed  in  step  3  is  term 
by  term  multiplication. 

The  tree  algorithm  routines  are  becoming  very  classical  and  therefore  we 
will  not  describe  them  either.  Note  however,  that  since  only  convolutions  are 
performed  in  our  algorithm,  we  only  use  the  discrete  Fourier  transform  of  the 
wavelet  coefficients  (and  not  the  corresponding  values)  at  every  scale,  that  re¬ 
duces  significantly  the  complexity. 

We  now  give  the  detailed  description  of  the  main  program  (ELLIP  )  and  of 
the  precalculus  program  (PRECAL)  in  pseudo  code. 

The  following  example  sketches  the  structure  our  programs. 

[OUTPUTS] =Pr ogram ( INPUTS ) 

#  Comments 

#  Body  of  Program: 

[0UTPUTS1]=  Subprogram! (INPUTS) 

INPUTS2  =  OUTPUTS! 

[OUTPUTS]  =  Subprogram2(INPUTS2) 

} 

Our  variable  descriptors  bears  some  resemblance  to  the  C  language  as  well 
as  to  the  MATLAB  conventions. 


III.4.1  Preliminary  computations 

The  symbols  and  ./  used  to  present  this  program  are  borrowed  from 

MATLAB  and  mean  respectively,  the  transposition,  the  matrix  product,  the 
term  by  term  product,  and  the  term  by  term  division.  We  also  use  the  following 
n  sub  sampling  operator  a  :  n  :  b  defined  as:  If  a  is  a  2D  array  of  size  2^  x  2^, 
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b(l  :  2^  :  2P,  1  :  2^  :  2^)  is  a  new  array  of  size  2^""”  x  2^’"”  given  by  a(i,  j)  = 
a{2-i,2Nj),{iJ)eIp. 

Program  P REGAL 

[QMFBIW,TAUTW]  =  PRECAL(p,pmax,QMFW,PHIW,SW,CW, SIGMA) 

#INPUT: 

#p  ">  index  of  the  approximation  space  Vp  in  which  the 

#  elliptic  problem  is  solve. 

#pmax  ->  index  of  the  approximation  space  Vpmax  in  which  the 

#  precomputation  of  TAUTILDE  is  done  (it  depends  on  the 

#  prescribed  precision) . 

#QMFW  ->  structure  containing  the  quadrature  mirror  filters  in 

#  one  dimension: 

#  QMFW.mO  ->  ID  array  containing  the  quadrature  mirror  filter 

#  coefficients  associated  to  the  scaling  functions; 

#  size(QMFW.mO)->2"p;  QMFW.mO(i)  =  m0(i/2>), 

#  i  belong  to  {0, .  . .  ,2"‘p-l} . 

#  QMFW.ml  ->  ID  array  containing  the  quadrature  mirror  filter 

#  coefficients  associeted  to  the  wavelet; 

#  size(QMFW.ml)->2"p; 

#  qMFW.ml(i)  =  ml(i/2^p),i  belong  to  {0, . . . ,2^p-l> . 
#PHIW  ->  ID  array;  si2e(PHIW)“>2"pmax;  where  pmax  is  given 

#  and  pmax>p;  PHIW(i)  =  the  value  of  the  Fourier  transform 

#  of  the  ID  scaling  function  at  the  point  i,  i  belong  to 

#  -Co, .  . .  ,2'"pmax--l}.  Used  to  compute  the  value  of  tautilde  on 

#  the  finer  grid. 

#SW  “>  2D  array  containing  the  sampling  of  the  function 

#  S  used  for  biorthogonal  filters; 

#  Size(SW)->(2^pmax  X  2^pmax) ;  S¥(i, j )=S(i/2^pmax, j/2"pmax) , 

#  (i,j)  belong  to  {0, . . . ,2"pmax-l>*2. 

# 

#CW  ->  2D  array  containing  the  sampling  of  the  function 

#  S(2w)/(2^s  S(w))  Size(CW)  ->  (2''p  X  2^p) ;  CW(i,j)  = 

#  S(2i/2''p,2j/2*p)/  S(i/2^p,j/2‘^p), 

#  (i,j)  belong  to  {0, . . .  ,2"P“1}‘'2. 

#SIGMA  ->  2D  array  containing  the  sampling  of  the  symbol 

#  of  the  operator  Size(SIGMA)  ->  (2*pmax  X  2"pmax); 

#  SIGMA(i,j)  =  sigma(i/2"pmax, j/2"p) , 

#  (i,j)  belong  to  {0, , . . ,2"pmax-l}"2 . 

# 

#0UTPUT: 

#QMFBIW  ->  structure  containing  the  biorthogonal  quadrature 

#  mirror  filters  related  to  the  tautilde  functions: 

#  QMFBIW.PTILDE0->  2D  tab  containing  the  biorthogonal 
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# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

#  Computation  of  the  filters  Ptilde 

# 

QMFBIW.PTILDEO  =  ((QMFW.mO)’  ♦  (QMFW.mO))  .*  CW; 

QMFBIW.PTILDEl  =  ((QMFW.mO)’  *  (QMFW.ml))  ./ 

(2‘'s  SW(l:2"(p-pmax) : 2*pmax, l:2''(p-pmax)  :2’"pmax); 
QMFBIW.PTILDE3  =  ((QMFW.ml)’  ♦  (QMFW.ml))  ./ 

(2''s  SW(l:2*(p"pmax)  :2"pmax,  l:2"(p-pmax)  :2"pmax) ; 

# 

#  Computation  of  the  point  value  filter  related  to  TAUTILDE 

# 

TAUTW  =  (  (2^(ps)  *  PHIW’  *  PHIW)  .*  SS)./  SIGMA; 

TAUTW  =  Periodize (TAUTW,  p) ; 

The  subroutine  Periodize  is  not  described  here,  but  it  is  a  straight  forward 
transcription  of  43. 

III.4.2  Main  Program 

Main  Program  EUip 

[UX]  =  Ellip(FX,QMFW,FIW,QMFBIW, TAUTW) 

#INPUT: 

#FX  ->  2D  array  containing  the  sampling  of  the  function 

#  f;  Size(FX)  ->  (2^p  X  2^p) ;  FX(i,j)  =  f  (i/2*p,  j/2'-p) , 

#  (i,j)  belong  to  {0, . . . ,2"p“l}*2. 

#FIW  “>  ID  tab  of  data  containing  the  interpolation  filter 

#  related  to  PHI_pO,  size(FIW)“>2"p; 

#QMFW  '>  see  PRECAL 

#QMFBIW  ->  see  PRECAL 
#TAUTW  ->  see  PRECAL 

# 


quadrature  mirror  filters  associated 
to  the  scaling  functions; 
size(QMFBIW.PTILDE0)->(2^p  X  2^p) : 

QMFBIW.PTILDE1*“>  2D  tab  containing  the  biorthogonal 
filters  associated  to  the  first  wav¬ 
elet  ;  size  (QMFBIW.PTILDEl  )~>  (2^p  X  2^p). 

QMFBIW.PTILDE2->  same  as  QMFW.PTILDEl  for  the  second 
wavelet  (not  computed  QMFBIW . PTILDE2  = 
QMFW.PTILDEl  transposed), 

QMFBIW. PTILDE3->  same  as  QMFW.PTILDEl  for  the  third 
wavelet . 
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# 

#OUTPUT; 

#UX  “>  2D  axray  containing  the  sampling  of  the  approxi- 

#  -mation  u_p;  size(UX)->(2"p  X  2"p),  UX(i,j)  = 

#  ii(i/2‘'p,j/2^p),  (i,j)  belong  to  {0, .  .  .  ,2^p-l}^2. 

# 

#TAUTW  ->  2D  array  containing  the  values  of  the  Fourier 

#  transform  of  the  scaling  function  TAUTILDE  at  level  p 

# 

#TEMP0RARY  DATA: 

#FW  ->  2D  tab  containing  the  fft  of  FX;  Size(FW)  -> 

#  (2^p  X  2^p); 

#CPW  ->  2D  tab  containing  the  fft  of  scaling  coefficient  of 

#  FX;  Size(CPW)  ->  (2^p  X  2*p); 

#DJW  ->  Structure  of  2D  array  containing  the  Fourier  transform 

#  of  the  wavelet  coefficients;  size(DJW)  ->  (2"p  X  2"p); 
#CTILDEPW 

#  “>  same  as  CPW  for  UX; 

#UW  ->  2D  tab  containing  the  fft  of  UX; 

# 

# 

#  step  0 

# 

[FW]  =  Fast  Fourier  Transform(FX) 

# 

#  step  1 

# 

[CPW]  =  FW.*FIW  (Term  by  term  multiplication) 

# 

#  step  2 

# 

[DJW]  =  2D  Tensorial  Tree  Algor ithm_D  (CPW,QMFW) 

# 

#  step  3 

# 

[DJW]  =  DJW.*(2‘'js)  (Term  by  term  multiplication) 

# 

#  step  4 

# 

[CTILDEPW]  =  2D  Non  Tensorial  Tree  Algorithm,.!  (DJW,qMFBI) 

# 

#  step  5 

# 
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[UW]  =  CTILDEPW .  ♦TAUTW  (Term  by  term  multiplication) 

# 

[UX]  =  Inverse  Fast  Fourier  Transform (UW) 

} 

III. 5  Storage  and  Complexity  Analysis 

As  the  computation  is  clearly  separated  into  precalculations  and  actual  imple¬ 
mentation  of  the  algorithm,  we  will  also  present  the  storage  and  complexity 
analyses  separating  the  two  parts.  One  should  remember  that  the  precalcu¬ 
lation  is  done  once  and  for  all  while,  as  it  will  be  the  case  in  section  IV,  the 
algorithm  can  be  applied  iteratively. 

We  will  not  discuss  the  complexity  related  to  one-dimensional  computations 
as  well  as  the  storage  connected  to  one-dimensional  arrays  since  both  can  be  ne¬ 
glected  in  our  bidimension al  implementation.  All  the  evaluations  are  performed 
for  N  ==  dim  Vp  =  2?^  x  2^. 

Storage 

Permanani  storage  (precalculations):  The  structures  QMFBIW  and  TAUTW  rep¬ 
resent  four  bidimensional  arrays  of  size  N, 

Temporary  storage  (actual  algorithm):  The  storage  related  to  bidimensional 
arrays  can  be  reduced  to  one  arrays  of  size  N . 

Finally,  the  total  memory  used  corresponds  to  five  arrays  of  size  A. 

Complexity  analysis 

Precalculus:  The  computation  of  the  four  arrays  in  the  structure  QMFBIW  is 
done  in  (7  X  iV. 

The  computation  of  is  performed  m  C  x  N  operations.  The  value  of 

C  depends  on  the  precision  of  the  calculation. 

Main  program: 

Fast  Fourier  Transform  and  Inverse  Fast  Fourier  Transform  in¬ 
volve  C  X  N  log{N)  operations. 

The  complexity  of  the  Term  by  term  multiplication  is  N. 

Tree  algorithm-D  and  Tree  algorithm-I  are  based  on  convolution  and 
decimation  operators.  These  procedures  involve  therefore  C  x  N  opera¬ 
tions. 

Therefore  the  total  complexity  is  0{Nlog{N)). 

In  the  following  section  we  use  these  programs  iteratively,to  solve  the  2D 
Burgers  equation  after  reducing  it  to  a  cascade  of  elliptic  problems.  We  would 
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like  to  emphasize  moreover,  that  our  approach  can  be  also  used  to  solve  equa¬ 
tions  involving  homogeneous  pseudo-difTerential  operators.  A  characteristic  ex¬ 
ample  is  \/-Au  =  /  with  periodic  boundary  conditions  on  [0, 1]^.  We  have 
L  =  and  therefore  (t(^)  =  Vir+W-  The  most  natural  choice  for  S  is 

S(0  =  2ysin^(^i/2)  -I-  sin^(6/2)  and  one  easily  checks  that  the  hypotheses  of 
theorem  II. 2  are  then  satisfied.  The  algorithms  presented  previously  can  be 
used  (see  Pj.  Ponenti  [16]). 


IV  Numerical  Application:  Resolution  of  the 
2D  Burgers  equations 


In  this  section,  we  will  use  the  periodized  Battle-Lemarie’s  multiresolution  anal¬ 
ysis  of  splines  of  order  m  (see  P.G.  Lemarie  [9]).  The  existence  of  collocation 
projectors  related  to  the  spline  breakpoints  requires  splines  of  even  order  and 
the  value  m  =  8  will  be  used  in  the  applications. 

As  described  in  J.  Liandrat  et  al.  [11],  any  parabolic  equation  of  the  type 


r  ^  +  Lou  +  G{u)  =  0 


u{0,t)  =  u{l,t) 
U  =  Uo  for  t  :=  0 


(44) 


(  0  <  t  <  Tmax^  0  <  a;  <  1 

where  Lo  is  a  differential  operator  of  even  order  wdth  positive  symbol  cro(u;),  and 
G  is  generally  a  nonlinear  function  of  u  and  its  derivatives,  can  be  numerically 
approached  using  a  classical  finite  difference  time  discretization  scheme  followed 
by  a  variational  approximation  of  the  resulting  elliptic  problems.  We  show  now 
that  the  approach  developed  in  the  previous  sections  can  be  used  efficiently  to 
provide  this  approximation. 

Following  J.  Liandrat  et  al.  [11]  we  first  introduce  a  segmentation  {tn}^=i 
of  [0,  Tmax]  (i.e.  a  sequence  such  that  0  =  f i  <  •••  <  =  Tmax) 

and  now  look  for  a  sequence  of  functions  of  the  x  variable  such 

that  is  an  approximation  of  u{x,tn)- 

With  Ain  -  fn+i  0  <  n  <  M,  and  considering  first  (44)  as  an  ordinary 
differential  equation  in  time,  a  standard  finite- difference  discretization  leads  to 
the  following  iterative  equation: 


£„«("+!)  =  f  (u(”), ...,  Af„, ....,  Giu(%  ...,  G(u("-^)))  (45) 

where  £„  is  a  step  forward  operator  that  together  with  F  is  determined  by 
the  choice  of  the  finite-difference  approximation  of  the  time-dependent  ordinary 
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differential  equation.  We  always  assume  that  this  approximation  is  at  least  semi 
implicit  for  the  linear  part  Lq  and  explicit  for  the  nonlinear  part.  Therefore 
3a  >  0  such  that  £„  =  (/  +  aAtnLo)  where  I  stands  for  the  identity  operator. 
By  hypothesis,  <to(uj)  >  0,Va;  and  then  Cn  has  always  a  symbol  bounded  from 
below  by  1. 

Hence,  assuming  that  and  =  0,  ...j}  are 

known,  the  resolution  of  (45)  falls  under  the  scope  of  paragraph  IIL3  and  the 
resolution  of  (44)  can  be  therefore  performed  iteratively. 

The  bidimensional  Burgers  equation  writes,  with  u  =  {ui,U2): 


^  4-  VuM  =  i/Au 
u{0,t)  =  u(l, t) 

u  =  uq  for  i  0 

^  0<t<  Tmax.O  <x<l. 


(46) 


Choosing  a  constant  step  segmentation  of  [O^Tmax]  (i.e.,  a  segmentation 
such  that  3At  such  that  VO  <  rz  <  M,tn  =  nAt),  an  implicit  Crank-Nicholson 
time  scheme  (second  order)  for  the  linear  term  (z/Au),  and  an  explicit  second 
order  Adams-Bashforth  scheme  for  the  nonlinear  term  (Vu.u)  we  get 


and  the  solution  can  be  written  as 


,(n+l)  3,  £^(n+l)  _  y(n)  ^ 


with 

{*("+!)=  (/-I/ 

To  fall  completely  under  the  scope  of  paragraph  III. 3,  one  should  be  able 
to  evaluate  the  point  values  of  the  nonlinear  term  of  (47).  We  used  the  sim- 
plest  method  available  that  consists,  as  classically  done  in  spectral  methods 
(C.  Canuto  et  al.  [3]),  to  “apply  the  nonlinear  operator  on  the  grid  points”. 
More  precisely,  the  approximation  of  we  used  is  PVCVu^^Ku^^^)  = 

Cy^{Cvv  X  (Vu^).Cv''7>(u^))  where  x  is  a  term  by  term  multiplication  of  finite 
sequences  and  Cy'p  is  the  collocation  projection  introduced  in  section  III. 2. 

Then,  for  each  time  step  nAt,  the  problem  clearly  belongs  to  the  cl2iss  of 
elliptic  problems  studied  in  section  (III)  with  L  =  I  -  i/^A  and  /  =  — 

The  iterative  form  of  the  equation  induces 
some  modifications  of  the  general  scheme  presented  in  III. 3  and  we  therefore 
provide  the  full  scheme  for  an  iteration  of  the  Burgers  approximation  scheme: 


2u(")  - 


(47) 
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0.  The  inputs  are  the  values  of  (u^’^^Jp)),  (^^(Jp))  and  (^^(-fp))- 

1.  F^'^\Jp)  -  2w(”)(Jp)  -  5(PF(Vu("\w(”^)(Jp)  is  computed  as  described 
above. 

2.  Fp"(x),  the  function  of  l^p^  interpolating  the  values  F"(.7p)  is  constructed 

as  fp”  =  <^vk^pk  > 

3.  Fp  is  then  decomposed  into  the  wavelet  subspaces  <  7  <  p  -  1  and 

as 

K=  E  (^'p",^a)^a+‘^oo 

AGAp' 

where  cqo  = 

4.  then  becomes 

4"+i)  =  Y1  2-^'''(^p">v>rK+coo 

AeAr' 

where  6^  =  {V^L-'^xp^Y . 

5.  up”'*’^^  is  then  expanded  in  terms  of  the  set  of  functions  {r^^,  k  e  Ip}  as 

=  E  ■ 

k&I^ 

We  also  get 

“  Eifee/p 

6.  From  the  point  values  of  Tj^{Jp),  ^tP{Jp),  and  ^Tp{Jp)  we  compute 
(y(n+i)(jp))  jtg  gradient  using  the  corresponding  point  value  filters. 

7.  Finally  using  the  values  of  (u„(Jp))  and  of  its  gradient,  we  get  the  values 

of  (m("+^^(Jp))  =  Jp)  -  ■«^”^(7p))  and  its  gradient. 

IV.  1  storage  and  Complexity  Analysis 

As  described  above,  the  numerical  code  implements  the  elliptic  solver  in  an  it¬ 
erative  process.  Since  the  time  step  At  is  constant,  the  characteristics  of  the 
elliptic  solver  does  not  depend  on  the  time  index  n.  Then,  the  solver  precalcu¬ 
lations  related  to  the  whole  parabolic  problem  are  the  same  as  the  ones  related 
to  the  elementary  elliptic  solver  (see  section  III. 5).  This  applies  to  the  storage 
and  to  the  complexity  as  well. 
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As  it  has  been  shown  in  the  previous  section  however,  extra  work,  not  con¬ 
nected  to  the  elliptic  solver  itself  but  to  the  computation  of  the  right  hand  side 
term  of  the  iterative  equation  (47)  is  required.  This  extra  work  is  related  to  the 
storage  of  the  fields  at  the  different  time  steps  involved  in  the  three  level  time 
step  Adams-Bashforth  Crank-Nicholson  scheme  and  to  the  point  value  evalu¬ 
ation  of  the  derivatives  involved  in  the  nonlinear  part.  Again,  it  can  be  split 
into  permanent  and  temporary  storage  as  well  as  in  precalculation  and  main 
program  extra  work. 

Storage  (in  addition  to  the  elliptic  solver  storage) 

Ptrmanani  storage  (precalculations):  One  extra  structure  containing  the 
point  values  must  be  stored  in  one  bidimensional  array  of  size  N;  the 

structure  containing  the  point  values  ■§^^p{Jp)  is  given  by  transposition  of  the 
previous  one. 

Temporary  storage  (actual  algorithm):  The  two  fields  and 

can  be  handled  using  three  arrays  of  size  N, 

Complexity  analysis  (in  addition  to  the  elliptic  solver  complexity) 

Precalculation:  The  computation  of  PV  a  is  performed  in  C  x  N  opera- 

dx  'p 

tions  where,  as  in  section  III.5  the  value  of  C  depends  on  the  precision  of  the 
calculation. 

Main  program:  The  addition  of  complexity  is  related  to  steps  (1),  (5),  (6) 
and  (7).  Since  these  steps  involve  convolutions  and  term  by  term  products,  the 
added  complexity  is  again  CNlogN 

Finally,  the  total  memory  used  is  7  arrays  of  size  N. 

The  total  complexity  is  0(iV’ log( A")). 

Obviously,  the  total  complexity  of  the  whole  resolution  is  M  times  the  com¬ 
plexity  of  one  time  step  resolution. 

IV. 2  Numerical  Results 

Test  case  on  an  0:2  translation  invariant  problem 

The  validation  of  our  code  has  been  performed  on  an  translation  invari¬ 
ant  problem  constructed  using  for  the  initial  condition  {uo^,uq^){xi,X2)  = 
(sin(27rxi),  0).  Indeed,  with  such  an  initial  condition,  the  solution  remains  x^ 
translation  invariant. 

For  an  easy  comparison  to  the  well  documented  paper  of  C.  Basdevant  et  al. 
[1]  we  used  1/  =  10“'^/7r. 

As  explained  in  C.  Basdevant  et  al.  [1],  the  pertinent  quantities  are 

dvi  du 

ms  =  sup(||^(®,t)||co  =  swp«e[o.Tm<2*]|^(0.5,<)| 
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and  tms  defined  such  that 


i^(0.5,^m5)|  = 
axx 

Table  48  exhibits  the  numerical  results  obtain  using  various  values  for  the 
time  step  At.  The  maximum  time  step  numerically  acceptable  was  At  =  0.0075. 
In  each  case,  the  values  of  ms  are  computed  by  interpolation  and  the  correspond¬ 
ing  values  of  tms  are  deduced.  The  comparison  with  the  expected  theoretical 
values  (first  column)  shows  that  our  method  competes  favorably  with  the  ma¬ 
jority  of  the  schemes  presented  in  C.  Basdevant  et  al.  [Ij.  A  complete  study  of 
the  time  step  size  dependence  of  the  results  connected  to  the  stability  analysis 
of  the  parabolic  algorithm  will  be  published  later. 


Exact 

A, 

0.0005 

0.001 

0.0025 

0.005 

-304.0103 

ms 

-304.6308 

-305.727 

-309.4354 

-316.5454 

0.255237 

^ms 

0.253 

0.252 

0.25 

0.245 

Test  case  on  a  first  diagonal  translation  invariant  problem 
Our  second  test  case  is  performed  on  a  first  diagonal  translation  invariant  prob¬ 
lem  constructed  using  (tioi , 2:2)  =  (sin(27r(a?i  -f-  0^2), sin(27r(a?i  +3:2)). 
Again,  the  solution  can  be  compared  to  the  reference  solution  of  C.  Basdevant 
et  al.  [1]  thanks  to  a  45"^  rotation  and  to  a  time  dilation  of  factor  2.  However, 
according  to  our  reference  axes,  it  is  obviously  a  fully  bidimensional  solution. 

Figures  1,2  and  3  show  the  isoline  values  of  the  numerical  approximations 
computed  with  6t  =  .001  at  i  =  0,  t  =  0.15  and  t  =  0.50.  The  first  diagonal 
translation  invariance  is  kept  and  we  obtain  the  values  ms  —  249.0528  and  tms  = 
0.123.  The  expected  theoretical  values  are  “304.0103  for  ms  and  0.1276185  for 
tms-  This  is  not  as  good  as  before  but  one  should  note  that  the  resolution  in 
the  direction  perpendicular  to  the  front  axis  is  now  half  the  one  in  our  previous 
calculations. 

Since  the  ultimate  application  of  all  this  work  is  the  development  of  adaptive 
algorithms  (i.e.  the  development  of  algorithms  handling  approximation  spaces 
of  reduced  dimension  adapted  to  the  solution  regularity  (see  for  instance  Pj. 
Ponenti  [16]), we  have  estimated,  at  various  times,  the  wavelet  basis  adapted 
to  the  approximation  and  defined  as  the  lowest  cardinal  m  =  8  spline  wavelet 
basis  preserving  the  norm  of  the  approximation  with  a  precision  of  lO""^. 
The  columns  of  table  (49)  show  for  each  scale  0  <  ;  <  7  the  number  of  wavelets 
selected  in  the  adapted  basis  related  to  the  approximated  solution  at  various 
times.  It  appears  that,  compared  to  the  full  basis  of  Vg  (last  column),  these 
bases  have  a  drastically  reduced  cardinal  (we  defined  the  rate  of  compression 
as  basis  ^  gradients  of  the  solution  fill  up  a  large 

domain  made  of  two  complete  lines  of  the  plane  (see  figures  1,  2  and  3) 
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Figure  1:  Initial  condition ,  ti  =  0. 


Figure  2:  Approximated  solution,  =  0.15. 
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Figure  3:  Approximated  solution,  ti  =  0.50. 
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Conclusion 

In  this  paper,  we  have  proposed  an  inversion  scheme  for  elliptic  problems 
based  on  biorthogonal  wavelets.  The  approximation  of  elliptic  problem  solutions 
is  constructed  and  leads  to  stable  and  fast  numerical  algorithms. 

Numerical  tests  related  to  the  approximation  of  the  parabolic  Burgers  equa¬ 
tions  transformed  into  a  cascade  of  elliptic  problems  are  provided. 

The  approximation  scheme  is  based  on  convolution  operators  and  can  there¬ 
fore  be  theoretically  used  in  the  framework  of  adapted  spaces  of  approximation. 
As  mentioned  however,  the  nice  tensorial  product  structure  that  enforces  nu¬ 
merical  efficiency  is  fragile  and  is  generally  destroyed  when  applying  the  scheme 
directly.  Other  approximations  for  the  step  forward  operator,  should  allow 
one  to  use  efficiently  this  approximation  in  a  general  context  of  adapted  multi¬ 
dimension  spaces  of  approximation. 


IN  umber  ol 
300  ti  =  500  wavelets  in  V7,- 
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